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A fluctuating non-ideal fluid at its critical point is simulated with the Lattice Boltzmann method. 
It is demonstrated that the method, employing a Ginzburg-Landau free energy functional, correctly 
reproduces the static critical behavior associated with the Ising universality class. A finite-size 
scaling analysis is applied to determine the critical exponents related to the order parameter, com- 
pressibility and specific heat. A particular focus is put on finite-size effects and issues related to the 
global conservation of the order-parameter. 



I. INTRODUCTION 

The theory of critical phase transitions has received 
considerable attention and seen remarkable progress in 
the past decades [D-Q- The importance of critical phe- 
nomena stems from the fact that systems whose micro- 
scopic behavior can be very diverse nevertheless share 
the same universal properties close to their critical point 
and thus belong to the same universality class. Uni- 
versality classes are defined only by a few characteris- 
tic properties, such as the number of components of the 
order-parameter, dimensionality of space, the couplings 
to other dynamical quantities in the system and the pres- 
ence of conservation laws. It should be noted that, uni- 
versality classes are different regarding static and dy- 
namic critical behavior: while, for instance, the uniax- 
ial ferromagnct and a binary fluid close to its demixing 
point both show Ising type static behavior, their critical 
dynamics is decisively different 

Most critical properties for standard bulk systems, 
such as critical exponents and amplitude ratios, are 
nowadays known with high precision due to the com- 
bined effort of experimental, theoretical and simulation 
approaches 0, @. Thus, in recent years, the focus 
has moved on to the study of critical phenomena in 
more complex situations, as, for instance, under non- 
equilibrium conditions [f|, at surfaces Q or in complex 
fluids [1, Q . Here, the hope is to utilize fluctuation in- 
duced effects, such as the the critical Casimir effect [Ioj |. 
for novel applications. Due to the increasing complexity 
of such systems, simulation approaches to critical dynam- 
ics in fluid systems thus become an indispensable tool. 

While dynamic critical phenomena of fluids have been 
extensively studied theoretically and by experiment d, d, 
ITll |. their simulation has onl y b een recently approached 
via Molecular Dynamics [12 - H5j . However, system sizes 
are rather limited and certain transport coefficients, such 
as the shear viscosity, are notoriously hard to determine 
with sufficient accuracy. Recently, the Lattice Boltz- 
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mann (LB) method - being a well-established and effi- 
cient solver of the Navier-Stokes equations - has been ex- 
tended to deal with thermal fluctuations in liquid-vapor 
systems [l|| as well as binary fluids [l7j]- Thus, the LB 
method appears to be a promising candidate for the sim- 
ulation of critical phenomena in simple and complex flu- 
ids. 

As a first step towards this aim, the current work 
presents simulation results on the static critical proper- 
ties of a non-ideal fluid obtained with the fluctuating LB 
model introduced in (l6| . The dynamical critical prop- 
erties will addressed in a separate paper [l8j . In the 
present model, the fluctuating hydrodynamic equations 
for the density and momentum of an isothermal, non- 
ideal fluid are solved via a Langevin approach. All simu- 
lations are performed in two dimensions, which - besides 
computational efficiency - has the advantage that crit- 
ical properties can be much easier assessed than in 3D 
since fluctuation effects are generally more pronounced in 
lower dimensions. Since the model is governed by a one- 
component Ginzburg-Landau </> 4 -free energy functional, 
the static critical properties are expected to be described 
by the 2D Ising universality class [19r{22| ■ This prediction 
is indeed borne out by the present LB simulations, which 
are also in line with previous Monte-Carlo investigations 
of the two-dimensional 4 -model [23Tj3l| . A crucial is- 
sue in a hydrodynamics based simulation approach is 
the global conservation of the order-parameter (here, the 
density) , which not only leads to a more pronounced criti- 
cal slowing down compared to the non-conserved case Q , 
but also complicates the finite-size scaling technique used 
to extract critical properties in a finite system |32l . HH . 
The aim of the present work is to provide a through as- 
sessment of the LB method in the critical fluctuation 
regime and demonstrate that, besides the above men- 
tioned complications, the method is capable of success- 
fully capturing critical fluctuations in fluids. At the same 
time, important issues that might be useful for further 
applications and extensions of the method shall be high- 
lighted. The paper is written in a self-contained manner 
and is hoped to provide also a researcher unacquainted 
with critical phenomena with sufficient background in- 
formation. 
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The outline of the paper is as follows. In section II, the 
critical properties of the Ginzburg-Landau theory are re- 
viewed. In particular, the effects of a finite system size 
and the global conservation of the order-parameter are 
discussed. Also, the order-parameter distribution as a 
fundamental quantity to extract information on critical 
as well as non-critical properties is introduced. Section 
III discusses the simulation method and contains a num- 
ber of remarks on the correct choice of simulation pa- 
rameters. In section IV, simulation results on the struc- 
ture factor and thermodynamic quantities are presented 
and compared to the theoretical predictions and previous 
works. 



II. THEORY 
A. Ginzburg-Landau model 

1. Introduction 

As the present simulation approach to critical phase 
transitions is based on fluctuating hydrodynamics (see 
section HTl]) , the fundamental quantity for our purpose is 
the density field p(r) of the fluid. From the density, an 
order parameter </>(r) can be defined as 



(r) 



P(r) - Po 
Po 



(1) 



where the reference density p is taken as the global aver- 
age, po = J drp(r)/V, with V being the system volume. 
The equilibrium behavior of the order-parameter is gov- 
erned by a Ginzburg-Landau free energy functional 



T[4>]= / dv -|V0| 2 + / o (0)-/# 



with /o being a Landau potential 



(2) 



(3) 



and h a possible external field, which, in the case of 
liquid-vapor criticality, can be taken as the difference 
p — p c between the actual and critical pressure. As usual 
k and u are strictly positive, while the coefficient r of the 
quadratic term can be either positive or negative, lead- 
ing either to a single minimum or a double- well form of 
the Landau free energy. Thermal fluctuations lead to a 
distribution of order parameter values according to the 
probability density 



z 



(4) 



which is a functional of the order-parameter field. In the 
above equation, k B is the Boltzmann constant and T is 
the temperature. Note that no dependence of the coef- 
ficient r in the Landau potential on the temperature T 



will be assumed. Rather, r is considered as an indepen- 
dent quantity representing the appropriate temperature 
measure in the context of Ginzburg-Landau models (see 
next section). The partition sum Z is given by 



Z = / V<f>e 



-T[<t>]/k B T 



(5) 



where J T>(/> denotes the integration over all possible re- 
alizations of the order parameter distribution. On a d- 
dimensional lattice of size S d with a total of N lattice 
points, the order parameter <j> is specified by its N values 
cf)(i) = 4>{ri) and the functional integral is regularized as 
/ T>(t> — > HfLi J d(j>i. Discrete e quiv alents for the deriva- 
tive operators can be found in |16j |. The partition sum 
([5]) allows to define a thermodynamic Helmholtz free en- 
ergy F and the corresponding density / in the usual way 
via 



F = fV = -k B T\ogZ. 



(6) 



From the free energy, eq. ^ , a global order-parameter 
M and susceptibility x can be formally defined as re- 
sponse functions with regard to the external field: 



M 



dh V 
d 2 f 1 



dh 2 k B TV 
V 



dr(<j)(r)) = (to) (7) 
drdr'[^(r)^(r'))-Wr))(^(r'))] 



k B T 



((to 2 ) - (to) 2 ) 



(8) 



where 



to = y J drift and the brackets denote aver- 
age with respect to the distribution P, i.e. (g (</))) = 
J T>4>g((f>)P{(f)} for an arbitrary function g of (f>. A fur- 
ther quantity of interest is the spatial correlation function 
of the order-parameter fluctuations 

C{r)5(v - r') = ((0(r) - <$)(0(r') - (0))) , (9) 

and its Fourier transform (structure factor) C(k). Note 
that translational invariance is assumed in the above 
equation. Related to the correlation function C is a 
non-local susceptibility \( r ) — C(r)/ksT, which can be 
defined analogously to eq. ([5} by considering the linear 
response to a spatially dependent external field. The def- 
inition of the specific heat, which quantifies the thermal 
response, requires some care, as a temperature change 
can be effected in several ways, depending on the param- 
eterization of the model. Here, the field theoretic con- 
vention 0, [2l| is followed and the specific heat is defined 
as the response with respect to a change of the coefficient 
r, 



ch 



a 2 / 

dr 2 



Ak B TV 
(0 2 (r))(^V))] 



drdr'[(0 2 (r)0 2 (r')) 



(10) 
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where E = J dr<j> 2 /2V represents the most singular part 
of the local energy T . 

It is often convenient to rewrite the Ginzburg-Landau 
free energy in terms of a minimal number of parameters. 
To this end, we note first that the temperature only ap- 
pears as an overall scale factor in the Boltzmann weight, 
cq. ((H) , and can thus be absorbed in the definition of the 
coupling constants. Second, the coefficient of the square- 
gradient term can be fixed to 1/2 by rescaling the order 
parameter field, <p = </>/\/ n/ksT. The reparameterized 
free energy functional then reads 



where 



T[<t>]/k B T 



I f |V0| 2 



uksT 



h(j) 
(11) 



(12) 



are the remaining independent coupling constants. Cor- 
respondingly, the Boltzmann weight in eq. @ becomes 
e -F{4>]_ <p ne function! (|Hj) is the usual starting point 
for field-theoretic studies of the Ginzburg-Landau model 
U HH HH. In the following, both parameterizations, 
eqs. © and (fTTf . of the model shall be used. 



2. Critical behavior 

In practical applications of Ginzburg-Landau theory, 
the free energy functional $2$ is obtained as a coarse- 
grained description of some microscopic degrees of free- 
dom, e.g. spins on lattice or molecules of a fluid (34l. l35l] . 
Typically, one defines the order-parameter as an aver- 
age of microscopic degrees of freedom over a certain 
coarse-graining length. In this regard, the free energy 
functional F can be considered as an effective Hamilto- 
nian from which a partition function and a correspond- 
ing Hclmholtz free energy can be obtained. Exact cal- 
culations starting from microscopic models often lead to 
more complicated effective Hamiltonians than described 
by the simple Landau potential of eq. ([3|). However, at 
least close to the upper critical dimension of the model 
(which is four in the present case) , standard renormaliza- 
tion group arguments show that all higher-order terms 
become irrelevant at the critical point and the simple </> 4 - 
free energy indeed describes the universal critical prop- 
erties of all systems with the same, Ising- type symme- 
try of the order-parameter [1 01, El III ■ Monte-Carlo 
simulations |23l — K31| as well as theoretical arguments in- 
voking conformal invariance [36T - l38j show that the two- 
dimensional </> 4 -model with a scalar order-parameter be- 
longs to the 2D Ising universality class. 

Neglecting thermal fluctuations and evaluating the 
partition sum ([5]) only along its saddle-point defines the 
mean-field approximation, for which the critical point oc- 
curs for f = r = 0. However, when the Landau poten- 
tial becomes very shallow, thermal fluctuations can eas- 
ily overcome the central potential barrier and therefore 



significantly contribute to the functional integral in (|5|), 
leading eventually to a breakdown of mean-field theory. 
The critical point of the full Ginzburg-Landau model now 
occurs at a slightly negative f , which, due to the non- 
linear interactions between the fluctuations, depends on 
the non- linear coupling u [!, H^] (see below) . The "dis- 
tance" to the critical point r c can be defined in terms of 
a reduced dimcnsionlcss temperature 



(13) 



where fixed n and T are assumed in the last equation. 
Above definition ensures that 8 > in the disordered 
phase (super-critical regime) and 8 < in the ordered 
phase (sub-critical regime). In mean-field theory, f c = 
r c = 0, thus definition (fT3"|) must be replaced by 8 = 
r/a, where a is a suitable constant in order to make 8 
dimensionless. 

Close the critical point, thermodynamic quantities typ- 
ically show a power-law dependence on the reduced tem- 
perature 8, with exponents that are identical for all sys- 
tems within the same universality class 0, |39[ . Two-scale 
factor universality implies that the singular dependence 
of the Hclmholtz free energy, eq. © , on the two relevant 
scaling variables temperature 8 and external field h is 
given by 



/ S i„ g (0,M Hfl| 2 - Q /±(VI0H 



(14) 



where f± is a universal scaling function (up to metrical 
factors) and a and 8 are critical exponents. From the 
above relation, the critical behavior of the order param- 
eter, susceptibility and specific heat follows as 



M ~ B{-8f (8 < 0) , 

x^r ± |0|-T, 

c H ~A ± \0\- a ^ A ± log(0), 



(15) 



where T±, B and A± are non-universal amplitudes (± 
refers to whether the critical point is approached from 
above or below). In the two dimensional Ginzburg- 
Landau model, the specific heat has a logarithmic diver- 
gence (which is conventionally indicated by an exponent 
a = 0) 1 . The correlation length £ diverges as 



£cx<r 



(16) 



while the correlation function at criticality assumes a 
power law, 



C crit (k) <x k~ 2 +\ 



(17) 



expected to be valid for k > l/£ |40H42j. The values of 
the critical exponents are collected in Table |U 



1 For other parameterizations of the temperature dependence of 
the model, there can also be a regular contribution to the specific 
heat. 
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Exponent 
Quantity 
Definition 


a f3 7 v rj 
specific heat order parameter susceptibility correlation length structure factor 

c H ~ <r Q m ~ e r3 x ~ d' 1 c ~ Q~ u c(k) ~ fc _2+ " 


Mean-field 
Ising 2D 


(disc.) 1/2 1 1/2 
(log.) 1/8 7/4 1 1/4 



TABLE I: Critical exponents for the mean-field and the 2D-Ising universality class. 9 is the reduced temperature. The specific 
heat is discontinuous in mean-field theory and logarithmically divergent in the 2D Ising case. 



From the reparamctcrized free energy, cq. fjl 1|) . we see 
that, in contrast to the Ising model, where only one cou- 
pling constant and thus a single critical point exists, the 
Ginzburg-Landau model entails a line r c (u) of critical 
points [j, |3{|. The universal critical properties of all 
points on the critical line are controlled by the renor- 
malization group fixed point, which is expected to be- 
long to the Ising universality class. The critical line of 
the Ginzburg-Landau model on a square lattice has been 
obtained in previous works via Monte-Carlo simulations 
[26U291 [3TI ] . Fig. [T] shows the corresponding phase dia- 
gram taken from [23]. Note that the critical value for f 
decreases with increasing interaction strength u. 

Two particular limits on the critical line deserve fur- 
ther remarks [43j]: In the "order-disorder limit", which 
is reached for r — > —00, u — > 00 with r/u = const., the 
potential has two minima separated by an infinitely high 
barrier and the lattice free energy functional becomes for- 
mally identical to the Ising Hamiltonian. On the other 
hand, the case r — > 0, u defines the so-called "dis- 
placive limit" , where the free energy functional is dom- 
inated by the gradient-term interaction and the central 
potential barrier is low compared to the thermal energy. 
This limit is particularly important in the case of struc- 
tural phase transitions [431 ]. In the infinite system limit, 
critical properties on all points of the critical line are uni- 
versal and of Ising type, except in the displacive limit, 
where Gaussian critical behavior is expected 0, H(| |43| . 
For finite systems, close proximity to the displacive limit 
leads to an undesired masking of Ising-type critical be- 
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FIG. 1: Phase diagram for the j 
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havior [li^. 

In the parametrization (|11[) of the Ginzburg-Landau 
model, dimensional analysis shows that the length di- 
mension of <fi is [<j)] = L 1_d / 2 , which immediately fixes 
the dimensions of the coupling constants as [f] = L~ 2 
and [u] = L d ~ 4 . Thus, in dimensions less than four, a 
dimensionless coupling constant can be defined as A = 
{t 2 /( 4 ~ d ) /f, which, in 2D becomes 



u u ksT 
f rn 



(18) 



since f and u have the same dimensions in this case. 
Fluctuation corrections to mean-field theory become sig- 
nificant for |A| > O(l). It should be remarked that re- 
lation ([L8|l is slightly misleading when used to charac- 
terize the critical region, as it seemingly implies that all 
critical points of the model lie on a straight line. As 
Fig. Q] shows, this is, however, only approximately cor- 
rect. In fact, a complete dimensional analysis must also 
take into account the existence of an additional length 
scale, e.g. a momentum cutoff A or a lattice constant 
needed to remove the infinities arising in the field theory. 
The presence of such a regulator then allows the critical 

= r c (uA d ~ 4: ) to become a non-trivial function of a di- 
mensional quantity like u. This fact is intimately related 
to the appearance of an anomalous dimension at the crit- 
ical point, leading to a scaling of the field as ^i-^/ 2 -?)/ 2 
under change of length. 

Below, some important analytical approximations to 
the Ginzburg-Landau model, that will be useful in an- 
alyzing the simulation results, shall be recapitulated 
briefly. 



3. Mean-field theory 

In the mean-field approximation, fluctuations around 
the order-parameter distribution <fi that globally min- 
imizes the Ginzburg-Landau functional are neglected 
This approximation underlies most non-ideal fluid 
LB models without thermal fluctuations and has been 
studied extensively in this context (see, e.g., 0, j4o|). 
The mean-field free energy is given by Fq = J-(<p) = 
Vfo((f>) and admits for two fundamental equilibrium so- 
lutions. One corresponds to a spatially uniform value of 
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the order parameter, given by 
- f 



{ ±yj-r/u 
and a mean-field susceptibility 



X 




(r > 0) 
(r < 0) 



(r > 0) 
(r < 0) . 



(19) 



(20) 



In addition, for r < there exists a solution describing 
an interface between the two free energy minima of the 
form 



with 



(z) = \J — r/utanh (z/w) 



w = (-2n/r) 1/2 



(21) 



(22) 



being the interface width. Note that, for the present def- 
inition of the interface profile, w is related to the mean- 
field correlation length £ [eq. (j2"o) ] by w = 2£. The sur- 
face tension (excess free energy per line in 2D) associated 
with the planar interface solution (|21[) is given by 



a= 3 



2nr s 



(23) 



For reference, the mean-field expression for the speed of 
sound is given by 



c s = P0- 



dp 2 



1 , 

— (r + 3uq 
Po 



— > (24) 
PoX 



where p is related to the global order-parameter 
eq. (P. 



4- Fluctuations 



by 




!(*) = 



FIG. 2: Expansion of the self-energy E in a self-consistent 
scheme up to 0(u 2 ). Thick lines represent the full correlation 
function C(q) . Dashed lines indicate amputated legs carrying 
the external momentum k. 



and c = —2 for r < 0. It is important to emphasize that, 
due to the assumption of translational invariancc, above 
expression for the structure factor for r < holds only 
in homogeneous states. In the above equation, 




(r > 0) 



(r < 0) 



(26) 



is the bare (mean-field) correlation length and x is the 
mean-field susceptibility, eq. (|20| , which is related to the 
correlation function as \ = h m k-»o Co(k) / ksT . 

In sec. Ill A 21 it was shown that the properties of 
the fluctuating Ginzburg-Landau model arc basically 
governed by the dimensionless coupling constant A of 
eq. (|18[) . It is informative to express this constant in 
terms of the physically more relevant parameters corre- 
lation length £, susceptibility Xi surface tension a and 
order-parameter <fi: 



X 



k B T X k B T(j) 2 



v 2 X 



(27) 



This shows that increasing the noise temperature T in the 
ordered state (small negative A) brings one always closer 
to the critical point, unless, for instance, the surface ten- 
sion or the density ratio ((j>) is increased accordingly. 



Thermal fluctuations around a uniform state can be 
systematically studied by splitting the order parameter 
into a uniform mean-field part and a spatially inhomoge- 
neous part </>(r) = <f> + 8<f)(r), where <j) = (<P) represents 
the average order parameter. Expanding T in the fluc- 
tuations 8(j> and treating the quartic anharmonicity as a 
perturbation allows to compute the correlation function 
as a scries of Gaussian averages, which can be conve- 
niently represented in terms of Fcynman diagrams. To 
zcroth order in the non-linear coupling u, one obtains 
the Ornstein-Zernikc expression for C (bare correlation 
function), 



C (k) 



k B T 



k R T 



1 



nk 2 



k £- 



k B Tx 

i + k 2 e 2 



(25) 



where the constant c accounts for the leading order con- 
tribution of the non-linear term to the correlation length 
in the phase-coexistence regime, being c = 1 for r > 



5. Perturbation theory 

In the critical region, the susceptibility is small and 
thus fluctuations have a strong effect on observable quan- 
tities, rendering the results of mean-field theory invalid. 
Not too close to the critical point, these effects can be 
described by perturbation theory, which shall be briefly 
summarized here 0, H^, H|| . The effect of the non- 
linear interactions between the fluctuation modes can be 
captured in terms of a self-energy S(k), which is defined 
by the resummed perturbation expansion ( "Dyson equa- 
tion" ) of the full correlation function C as d |2j [H |46| 



C7(k) 



1 



C - i (k)+S(k) 



(28) 



£ is given by the sum of all two-point one-particle irre- 
ducible diagrams. In the symmetric phase, the diagram- 
matic expansion of the self-energy in the self-consistent 
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scheme up to second order in the coupling u is shown in 
Fig. [2J where the solid lines represent the full correlation 
function C. The one- and two-loop-contributions to the 
self-energy are given by 0, [U [22|, |46{ 



E (i) 



3u 



^o(q) 



£< 2 )(k) 



(27T) 



d d qi d d q 2 



(29) 



fc s T y J (2n) d {2ir) d 
Go(qi)C (q2)Co(k- qi -q 2 ). (30) 



The notation J indicates that the integral has to be cut 
off at a wavenumber A. In the present case, the cut-off is 
provided by the lattice constant and above integrals have 
to be understood as sums, 



d d q 

(2n) d 



L d Yj 



(31) 



The sum runs over all permissible wavevectors on the 
lattice except the zero-mode, which must be excluded 
owing to the global conservation of the order-parameter. 
Note that T,^ is independent of the external momentum 
k. 

In a self-consistent treatment, the full correlation func- 
tion C is taken to be of the same form as Co but with 
renormalized parameters r' , k', i.e. 



C(k) 



K'k 2 + 0(k 2 



(32) 



Since £ itself depends on the renormalized r' and 

eq. 



represents a system of two coupled integral equa- 
tions for the determination of r' and k' from the bare 
parameters r and K. The wavevector-independent part 
of the self-energy, S(0) obviously renormalizes the sus- 
ceptibility parameter r, 



fc B TE(0) 



(33) 



The wavenumber-dependent part of S, which is of two- 
loop order, renormalizes the square-gradient parameter k 
and ultimately gives rise to a non-zero anomalous dimen- 
sion rj at the critical point. At criticality, S(k) scales as 
k 2 ~ n . Analogously, the fluctuation contributions to the 
coupling constant u can be determined from the vertex- 
corrections to the four-point correlation function, which 
are also at least of two-loop order. Taken together, one 
obtains a system of three coupled integral equations for 
r' , k' and v! in dependence of the bare parameters. For 
the present purposes, however, it is sufficient to focus 
only on the dominant effect, which resides in the renor- 
malization of r. Of course, this limits the predictions to 
a region not too close to the critical point. Using cq. (f2"9")l 



and (|30|) , the self-consistency equation for the renormal- 
ized temperature parameter r' follows as 

r = r + 'iuk B T [ 777-TT-r-- — 2 ~ 6u 2 (fcsr) 2 
J {2n) d r + nq 2 

d d qi d d q 2 1 1 1 

(27r) d (2n) d r' + Kqf r' + nql r' + «(qi + q 2 ) 2 

(34) 

which can easily be solved numerically. In practice, 
eq. (|34|) is used to find, for a given r employed in a simula- 
tion, the corresponding value of r', which will then allow 
one to compute the physical (renormalized) susceptibil- 
ity X — l/ r ' an( i correlation length £ = (ft/r') 1 / 2 . This 
will give sufficiently accurate predictions in the crossover 
regime from mean-field to the critical region to be com- 
pared to simulation results. 

Below the critical point, the order-parameter acquires 
a non-zero expectation value (0) , which, at the mean- field 
level, is given by eq. (|19[) . Fluctuation corrections, how- 
ever, induce a shift of the mean-field expectation towards 
smaller values. This effect can be isolated by splitting the 
order-parameter as 

^(v)=v + a(r), (35) 

where v = ((f)) is enforced by requiring a vanishing ex- 
pectation value of the fluctuation [U |48| 



M = . 



(36) 



Inserting eq. ([33]) into the free energy functional © leads 
to (up to an unimportant constant) 



= / dr 



f|V.| 



2 V 



3uv 2 )a 2 



+ (rv + uv 3 )a + uva 3 H — ua 4 
4 



(37) 



The last three terms can be considered as a perturbation 
around the Gaussian part given by the terms quadratic 
in a [13, S3 2 . To first non-trivial order, eq. (|3"o]) is rep- 
resented in diagrammatic form by Fig. [3^, and follows as 



= ((f)) = rv + uv 3 + 3uv 



d d c 



(2-K) d Kq 2 + (r + 3uv 2 ) ' 
(38) 

which defines an implicit equation to be solved for the 
true v. Note that the first two terms lead to the mean- 
field result for v, vmf = (—r/u) 1 / 2 , while the last term 
gives the first-order fluctuation correction. The corre- 
lation function of the shifted field a is obtained from 
cq. (f37"l) as 

k B T 



a(k) 



Kk 2 



fc B TE ff (k) ' 



(39) 



2 To the order of perturbation expansion that will be considered 
here, a distinction between the bare and renormalized r in l|37( 
is not necessary 
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{*}= x + 

(a) 





FIG. 3: Perturbation theory in the broken phase: (a) Contri- 
butions to the expectation value of <f> to 0(u). (b) Contribu- 
tions to the self-energy £(fc) to 0{u). In the broken phase, a 
new 4> 3 interaction with a coupling uv appears. Note that, as 
v ~ it~ 1//2 , the first contribution to £(fc) is in fact of 0(u). 

where r a = r + 3uv 2 represents the inverse bare suscep- 
tibility of a and the leading-order self-energy corrections 
are given by the diagrams in Fig. [3}d, amounting to [46l — 

m 

+ 3«/ 7?%^^- (40) 
7 (27r) d K q 2 + r a 

From eq. (|39|l one obtains the true, renormalizcd suscep- 
tibility x' = with 

r' a = r a + fc B T£ CT (0). (41) 

Analogously to the situation in the symmetric state one 
could increase the accuracy of the perturbation expan- 
sion by replacing all appearances of r a in the self-energy 
S CT (0) by r l a% thereby taking implicitly into account the 
fluctuation corrections to the correlation function given 
by the diagrams in Fig. [3)3 to all orders. However, to 
the order of perturbation theory set up in eq. (|1D|). the 
difference between the two expressions is negligible. 

B. Finite-size effects 

On approaching the critical point in an infinite sys- 
tem, various intensive thermodynamic quantities display 
power law divergences (see Table IJ). In a finite system, 
any quantity must necessarily stay finite and the critical 
divergences appear rounded [jj, l50( . Typically, deviations 
from the true critical behavior set in once the correlation 
length £ ~ 0~ u of the infinite system formally exceeds the 
system size S. In this case, one enters the so-called finite- 
size scaling (FSS) regime, where the physical correlation 
length scales with the system size S. Standard finite-size 
scaling theory [l], , which shall be summarized here, 
asserts that in this regime the power-law dependence on 
the temperature of a thermodynamic observable O in the 
infinite system, O ~ Q~ x , essentially transfers to a power- 
law dependence on the system size 

O ~ S x / v ~go{S/0 ~ S x '"g o (S 1 ' v 0) , (42) 



where go and go are universal scaling functions 3 . To 
ensure that the correct asymptotic limit for the infinite 
system is reached, one must have go{ z ) ~ ^4-e>±M _:c a s 
z — > ±oo [where Ao± denotes the corresponding ampli- 
tude, cf. eq. (fl"5j) ]. while go(z) must be regular for z — > 0. 
Note that, if O represents the order-parameter, only the 
limit z — > — oo is relevant, since the order-parameter is 
zero in the symmetric phase. For the specific heat in 2D, 
FSS theory predicts that 

c H ~]og(S)gc(S 1 ' v 6). (43) 

In a finite system, the maximum r Cj o(L) of an observable 
O defines an apparent critical point, which is typically 
found at a slightly different temperature r than the true 
critical point r c of the infinite system. The latter can be 
inferred by extrapolating the apparent critical point to 
the limit S — > oo using r c ,o{S) = r c (l + aoS~^' 1 '), with 
a constant ao For the infinite system, all apparent 
critical points must merge. 

In the present simulation approach, the density is glob- 
ally conserved, implying that the global order parameter 
ms = Yli 'Pi an d derived quantities, such as the suscep- 
tibility (compressibility), xs = fefr(( TO s) ~ ( m s) 2 ) > arc 
trivially zero. Therefore, a standard FSS study based 
on the total system size S, as in eq. (|4"2"j) . is not possi- 
ble in this case. Instead, ideas originally developed for 
grand-canonical simulations of the Ising-model [HI, HIJ - 
which have been later successfully applied to canonical- 
ensemble simulations of lattice gas models and off-lattice 
fluids [lil, |33| - shall be followed here. These methods 
essentially consist of dividing the total system of length 
S into subsystems (blocks) of smaller length L = S/2 1 , 
for integer i 4 , and computing the quantities of inter- 
est in these subsystems. In particular, a coarse-grained 
order-parameter can be defined as 

ieb 

where i runs only over the lattice nodes that lie in the 
given subsystem b. Note that for L = 1, each block 
corresponds only to a single lattice site and thus 
becomes identical to the field variable 4>{r^) at that site. 

While the coarse-grained order-parameter m£ exhibits 
fluctuations, its average {mi) only gives information on 
the global asymmetry between amounts of the two phases 
that are present; in particular, it can still be zero even in 
the phase coexistence regime. Thus, instead the quantity 

M L = (\m L \) = -^<|][>|) (45) 

ieb 



3 A universal scaling function can depend on its argument as 
f(a.Qz) where ao is a non-universal constant [l[ 

4 In principle, one could also allow arbitrary integer divisions of 
the system size instead of powers of 2. 
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shall be considered as the appropriate block order- 
parameter for all temperatures, in agreement with the 
convention employed in the Monte-Carlo method [531 ]. 
Note that, in a finite system, Ml will be non- vanishing 
even in the disordered phase, but will approach zero in 
the thermodynamic limit L — > oo. In eq. (|45[) . the aver- 
age is performed over all blocks b with the same size L 
and over time. 

To define a susceptibility based on the average order- 
parameter, the disordered and ordered regimes have to 
be considered separately [Hj|. In the disordered phase, 
the standard definition 



XL 



k R T 



« 



> 0) (46) 



is employed. In the ordered phase, a slightly modi- 
fied definition has to be used to ensure that xl only 
measures fluctuations around the equilibrium order- 
parameter value: 



XL 



((m|) - {\m L \? 



(6 < 0) 



(47) 



For deep quenches into the ordered regime, pronounced 
interfacial effects prohibit a direct application of defini- 
tion (|47p. In these cases, it is found here that the interfa- 
cial contributions have to be explicitly removed from the 
underlying order parameter distribution (see below) in 
order to obtain a reliable estimate for the susceptibility. 
In the non-critical regime, where the correlation length is 
<C L, above subbox susceptibility disagrees with the true 
susceptibility (as obtained, for instance, from the corre- 
lation function) by a boundary correction proportional 
~ £/L [H, [33[ • This can be seen by rewriting relation 
(Hi as 



XL 



i 



k B TL d 



(48) 



and noting that, if i is close to the subsystem boundary, 
possible correlations between <f>i and the order-parameter 
outside the particular subsystem (but within a range of 
£) arc not taken into account in the sum over j. Thus, 
the subsystem susceptibility xl deviates from the true 
susceptibility by a surface-to- volume ratio ~ 1/L. In 
complete analogy to the susceptibility, a coarse-grained 
specific heat can be defined as 



ch.l 



L d 



{{El) {E L ?) 



(49) 



where 



is the average of the energy-like parameter in a subbox. 



When applied to the block obscrvablcs Ol, the original 
FSS ansatz, eq. P^j). must be extended by an additional 
scaling variable L/S [33j] . which is necessitated by the 
fact that for L = S order-parameter fluctuations are ab- 
sent and that, consequently, corrections to scaling are 
expected to depend on the ratio L/S. Thus we can write 



Or 



L x l v go{L L l v Lj S) , 



(50) 



and similarly for the specific heat. Strictly, the above FSS 
ansatz is expected to be valid for 9 — > and 0<L<S 
with L — > oo, while, outside this range, the influence of 
further corrections to scaling (for example, induced by 
the presence of irrelevant scaling fields) will become no- 
ticeable [5(1 H3- However, in the case of Monte Carlo 
simulations, it is often found that the simple FSS rela- 
tion (|50p works surprisingly well already for rather small 
lattice sizes [52|. Relations ([50)1 will therefore be used in 
the following for analyzing our simulation results. 

The utility of the above FSS relations is based on the 
fact that they allow for a determination of universal crit- 
ical parameters, such as exponents and amplitude ratios, 
provided that the location of the critical point is known. 
If the critical point is not known, one might still obtain 
reasonable estimates for the exponents and the critical 
temperature by trying different values until a good scal- 
ing of all data points is achieved. Relation eq. (|5T)]) is par- 
ticularly useful in a case where the exponents are already 
known and instead the critical point has to be located. 
In this case, one plots OL~~ x l v for different L versus the 
coupling u for a fixed value of r. By eq. (|50[) we have 
(neglecting the dependence on L/S) 



(51) 



and thus, all curves cross through the same point when 
the critical coupling u c is passed [U [55| . 



C. Order- parameter distribution 

Quantities such as the coarse-grained order parame- 
ter or susceptibility can be generally defined from the 
moments of an underlying order-parameter distribution 
function Pl corresponding to a given subsystem size L 

MM, 



(m L ) = / m k P L {m)dm . 



(52) 



The distribution function is particularly useful in the crit- 
ical regime, where the order-parameter fluctuations have 
a pronounced non-Gaussian character that can not be 
fully captured by the low-order moments of Pl alone. In 
the phase-coexistence regime, the distribution moreover 
contains crucial information on surface tension [52l . [56j 
and phase-equilibria [57j . Through its intimate connec- 
tion to a coarse-grained free energy, it is also of funda- 
mental importance in the description of nucleation and 
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spinodal decomposition processes [57H6G |. Practically, 
Pl is obtained in a simulation by creating a histogram of 

the coarse-grained order-parameter from all subsys- 
tems. The above order-parameter distribution is there- 
fore a coarse-grained quantity, and should thus not be 
confused with P[<f>] of eq. (j4j, which is a functional of 
the order-parameter field <fii, i.e. P[<fi] = P[{<j>i}] depends 
on all values of cj> on the lattice. Formally, the coarse- 
grained distribution P L can be defined as a constrained 
average over all order-parameter fluctuations compatible 
with a given value of the average order-parameter tul in 
a cell of size L, 



cial effects, as argued in 33J. The cumulant analysis will 



P L {m) 



i 

T d 



-r[<j>]/k B T 



(53) 

where Z is a normalization factor. For instance, for 
L = 1 the distribution Pl is obtained by integrating 
the probability functional P[{0.;}] over all <\>i except one, 
Pi(^)=IUi tj dct>iP[{<f>i}]- 

The connection of the constrained distribution Pl (to) 
of cq. (|53l) to thermodynamics can be made explicit by 
definin g a constrained Helmholtz free energy Fl (to) via 

F L (m) = -k B T\og{P L (m)Z) = -k B T log P L (m) + F , 

(54) 

where relation ((6]) has been used. The thermodynamic 
Helmholtz free energy F, eq. ©, is obtained from Fl by 

/+oo 
dmexp[-F L (m)/k B T}. (55) 
-oo 

It must be emphasized that Fl is in general different from 
the bare Landau potential fo [eq. ([5])]. as the former is 
derived from an integral over the full probability func- 
tional and therefore includes (due to the gradient term) 
the effects of interactions between the fluctuations. Ex- 
ceptions are the limit T — > as well as F\ in the Ising 
limit, since in these cases the interaction between dif- 
ferent cells can be neglected. In the large volume limit, 
the constrained free energy can be shown to be equiv- 
alent to the effective potential often employed in field 
theory [64|. Above the critical point, the constrained 
free energy is a direct measure for fluctuations of the 
coarse-grained order-parameter wll around the equilib- 
rium state. Below the critical point this picture breaks 
down, as Pl will then not only receive contributions from 
homogeneous fluctuations but also from the presence of 
two-phase states. Nevertheless, Fl is often found to be 
well approximated by a simple Landau form, i.e. a low- 
order polynomial in m [62l [63l l65T - [67j . An often invoked 
alternative characterization of the order-parameter dis- 
tribution that avoids fit ting a potential is based on its 
higher-order cumulants [52|. However, it is found that, 
in the present case, the standard cumulant ratio U3 does 
not appear to have a well-defined L-indcpcndent limit at 
the critical point. This is most likely caused by intcrfa- 



therefore not be pursued further in the present work. 

The shape of the coarse-grained distributions Pl and 
hence the free energies Fl can be anticipated based on 
simple physical arguments [HI, |(35[ . Far above the criti- 
cal point, non-linear effects are small and thus the order- 
parameter distribution is expected to be well approxi- 
mated by a Gaussian centered around the average order- 
parameter value (m) = (52j . 



P L {m) 



1 



(27r<m£))V2 

L d / 2 
{2^k B T XL ) 1 l 2 



exp 



2(to|; 



exp 



It d 



m z L 



2k B T XL 



(56) 



where relation (|46|) for the variance (mi) 2 has been used. 
The width of each Gaussian decreases with larger coarse- 
graining length L as more and more fluctuations are 
averaged out and the distribution approaches the high- 
temperature fixed-point. 

Below the critical point, the shape of Pl depends dis- 
tinctly on the size of the coarse-graining length L in com- 
parison to the correlation length £. For the case of a 
phase transition at the critical density (that is exclu- 
sively considered here), the global conservation of the 
order-parameter requires that below criticality equal vol- 
umes of liquid and vapor are present in the simulation 
box. For L < £, a given subbox will typically cover ei- 
ther liquid or vapor and thus Pl is expected to show 
two approximately Gaussian peaks centered around the 
spontaneous values of the order-parameter ±Ml [which, 
to a first approximation are given by eq. (|19p]. 



1 



P L {m) = - 



L d/2 



2 (2nk B T XL ) 1 / 2 



exp 



exp 



(m 



(to - M L ) 2 L d 

2k B T XL 
M L ) 2 L d 



2k B T XL 



(57) 



In general, the region between the peaks of the distribu- 
tion, — I Ml I < m < I Ml I, represents not only the proba- 
bility of homogeneous order-parameter fluctuations, but 
also of the occurrence two-phase configurations in a sub- 
system. For L < £, the box cannot cover complete 
phase-separated states, and thus for this case the height 
of -Pl(O) is essentially a measure for the probability of 
homogeneous fluctuations. In general, however, homo- 
geneous fluctuations are exponentially suppressed by the 
volume L d [eq. ([5"7])] , in contrast to heterogenous fluctu- 
ations, whose free energy cost is just proportional to the 
area of the interface, Thus, for subsystems with 

L>(, homogeneous fluctuations give a completely neg- 
ligible contribution to the central region of Pl, and one 
can estimate the probability for a hetcrophase fluctuation 
as 



Pl(0) — constL 21 exp - 



2L 



d-l, 



k B T 



(58) 
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with an empirical exponent x that is typically found to 
be close to zero [5(1 [68|]. In the ensemble average, liquid 
and vapor phases will occur equally often and with any 
proportion in each subbox. Thus, one expects that for 
large L the central region of Pj, will become flat and, in 
the limit L — > oo, where interfacial contributions become 
negligible, eventually attain the same level as the peaks. 
This is also expected based on the notion of a coarse- 
grained free energy Fl , as for L — > oo the thermodynamic 
limit is reached where the true free energy must be 
convex by thermodynamic stability (and implied by the 
Maxwell construction) [6^, [7(| . 

In the vicinity of the critical point, i.e. in the finite-size 
scaling region characterized by £ > L, the distribution 
function becomes markedly non-Gaussian and one can 
make a scaling ansatz as 

P L (m) = L^ v P{mL / u , L X I V B, L/S) (59) 
where P is a universal scaling function 

[Em, mis As 

can be easily checked, above relation reduces to the cor- 
responding FSS relations for the moments, cq. (|50j) . if ad- 
ditionally use of the hyperscaling relation dv = 7 + 2/3 is 
made. It must be emphasized that the scaling function P 
is universal only for sufficiently large L, where it embod- 
ies the collective features of the critical phase transition. 
For small L, in fact different shapes for Pl at criticality 
are possible, depending on the location on the critical 
line. Towards the displacive limit, the barrier between 
the minima of the local potential is low, leading to a 
nearly Gaussian shape of Pl , in contrast to a pronounced 
double-peak structure in the Ising limit, where Pi closely 
reflects the on-site potential /o LLL( . The large-scale prop- 
erties of the order-parameter distribution (in 2D and 3D) 
have been extensively studied in previous works via field 
theoretic approaches [6l|, [7lT473l and Monte- Carlo simu- 
lations of Ising- like systems [U [6^, [H, [H, [66|, [7!| . 



III. SIMULATION METHOD 
A. Model 

In the present work, the equilibrium behavior of the 
Ginzburg-Landau model is simulated on a square lattice 
via a fluctuating hydrodynamics approach. Specifically, 
the Langevin extension of the non-ideal fluid LB model of 
Swift et al. [zl, [z3, introduced in [l6|, is employed. The 
LB equation (LBE) can be understood as a discretization 
of the continuum Boltzmann equation and contains as a 
subset the Navier-Stokes equations in the limit of long 
time and length scales [78|. The LBE describes the evo- 
lution of a set of distribution functions fi(r) = /(r,Ci) 
on a lattice streaming along a finite number of possible 
velocity directions Cj linking the nodes. Here, simula- 
tions arc performed on a D2Q9 lattice, i.e. the space 
dimension is two and i = 1, ... ,9. Employing a simple 
BGK-approximation to the collision operator, the present 



LB model is defined by the evolution equation 

fi(r + Ci At,t + At) 

= /j(r,t) - — [fi{T,t) - /<>,t)] +&M) . (60) 

T 

where t is the time, At is the time step, r is a relax- 
ation time, / 4 eq is the equilibrium distribution and is a 
random force term, to be specified below. The relevant 
observable quantities are given by the low-order moments 
of the distribution. In particular, we have for the density 
p and the fluid velocity u, 

p = I>=E/rs ^=£/<* = ( 61 ) 

i i i 

In the model of Swift et al. [76|, [77}, equilibrium ther- 
modynamics as embodied by the Ginzburg-Landau free 
energy functional is implemented by requiring the sec- 
ond moment of the equilibrium distribution to recover a 
thermodynamic pressure tensor P, 

£ CiaCipfP = P a /3 + pU a U/3 
i 

+ v{u a d(jp + upd a p + Uydjpdap) . (62) 

The term proportional to the kinematic viscosity v = 
v/p = ( r — l/2)/3 is introduced in the above equation to 
improve Galilean invariance [79|. The pressure tensor P 
is given by 

Pa/3 = (po - npV 2 p - ^|V,o| 2 ) 8 a p + K(V Q p)(V ( gp) , 

(63) 

where po = p-^d^fo — fo is the thermodynamic pres- 
sure. The pressure tensor satisfies the relation V ■ P — 
pW(5F/5<p) and can be obtained from the free-energy 
functional ([2]) , for instance, via the Noether theorem, the 
principle of least action or from the requirement of hy- 
drostatic equilibrium [80l - [83| . Physically, it accounts for 
the energetic balance between changes in fluid structure 
due to advection and surface-tension [§3] and thus en- 
sures that the equilibrium order-parameter distribution, 
eq. (U), remains unchanged by the flow [35[. An explicit 
expression for the modified-equilibrium distribution / eq 
on a D2Q9 lattice can be found in [Hj]. 

To complete the description of the employed LB model, 
the properties of the noise variables £j have to be spec- 
ified. In order to properly account for mass and mo- 
mentum conservation, the moment representation of the 
LBE is invoked [HI, H3| ■ This representation is defined 
by a set of basis vectors T ai (a = 1, ... ,9) that admit 
the distribution function to be expanded in terms of 
its moments m a as 

711 ■ 

fi{r,t)=T ai j±m a (r,t). (64) 

Here, Wi are a set of lattice weights and N a are the 
squared lengths of each basis vector T a . In the present 
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work, T a i is taken to be the basis set given in [T(j. The 
first three moments m a then follow as p,pu x , and pu y , 
while the higher moments cover the stresses (a = 4, 5, 6) 
and so-called ghost-modes (a = 7, 8, 9). The moment rep- 
resentation of the noise is now simply given by £ a = T a j£j, 
and only the result for £ a shall be stated below, which has 
a much simpler expression than £j . As was shown in [l6j , 
due to the use of a modified-equilibrium distribution that 
incorporates the non-ideal gas thermodynamics (|63l) , the 
noise obtained from the exact fluctuation-dissipation the- 
orem is wavenumber-dependent. This is clearly undesir- 
able, as such a form of noise is not directly applicable 
to spatially inhomogencous situations involving phase- 
separation. However, as the offending terms in the noise 
covariance are proportional to the square-gradient pa- 
rameter k, it is possible, by reducing k appropriately, 



to employ spatially uncorrclatcd noise while still main- 
taining satisfactory equilibration. Additionally, spatially 
uncorrelated noise has the advantage of being straight- 
forward to implement and computationally cheap. 

The noise is thus taken to be a Gaussian random vari- 
able without explicit correlations in space or time 

(£a(r, t)i b (r', t')) = E ab {r)S r , Tl S t ,p . (65) 

However, in order to properly account for spatially inho- 
mogencous fluid properties, occurring, for instance, in the 
ordered regime, the covariance S is still allowed to depend 
locally on position |88| . The final expression for the noise 
covariance, taking into account above-mentioned modifi- 
cation and neglecting terms proportional to the square- 
gradient parameter k, is obtained as [l6j 



_ 3p(v)k B T 1 / 1 
" ( ' ~ AVAt t 1 



4[2-3cg(r)] 



4/9 



12[c*(r)-|] 



(66) 



1/9 



2/3 



12 [c»(r) - §] 



2/3 



16 [§-§ c 2(r)]/ 



Here, AV is the volume of a lattice cell, which, together 
with the time step At is equal to one in lattice units. The 
quantity fc^T fixes the fluctuation amplitude and is es- 
sentially a free parameter of the LB model, only subject 
to the low-Mach number constraint of LB. Above noise 
covariance ensures that that fluctuations of the fluid ve- 
locity obeys locally the equipartition theorem of statisti- 
cal mechanics 

(u a (r)u^(r')) = —r--5 al3 5 r , r , . (67) 

In a simulation, Gaussian noise with a non-diagonal co- 
variance matrix can be created via a Cholesky-transform 
(see 0). 

As a consequence of the LB dynamics of eq. (f60|) . at 
large length and time scales the density and momentum 
obey a continuity equation 

d t p + V ■ (pu) = . (68) 

and a momentum-conservation (Navier-Stokcs) equation 
for a non-ideal fluid, 

dt(pu) + V ■ (pun) = -V-P + V- ct + V- R, (69) 



where 

<J a f3 = V (d a up + df)U a — -^d~t u ~^J + (dyiiy (70) 
is the viscous stress tensor, 




is the shear viscosity (which should not be confused with 
the anomalous dimension index) and 

C = f (r - (2 - 3c^) (72) 

is the bulk (or volume) viscosity. Noteworthy, in the 
BGK-approximation to the Boltzmann equation, the 
(bare) viscosities depend on the local density [89| . In the 
critical regime, however, it always possible to choose the 
parameters in the Landau free energy such that the mag- 
nitude of the density fluctuations is small compared to 
the background density, Sp/p 1 (cf. Fig. [To]) . If the vis- 
cosities are approximated as constants, and furthermore 
the non-linear advection term (which is not relevant at 
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criticality [4|) is neglected, the Navier-Stokes equations 
simplifies to 



d t (pu) = -V-P+7yV 2 u+(C + [1 - 2/d]r)) VV-u+V-R , 



(73) 



Below the critical point, where appreciable density vari- 
ations are present, the full cq. (|69|) must be used. Note 
that expression (JT2J) for the bulk viscosity differs from 
the standard LB expression by a factor (2 — 3c 2 ), which 



J 



(i? a/3 (r,*)i? 7(5 (r',t')) = 2k B T 



is an artifact of the modified equilibrium distribution of 
the present LB model filll77| 5 . 

The random stress tensor R imparts thermal noise on 
the fluid momentum, which is then transferred to the 
order parameter sector, leading - in equilibrium - to 
thermal fluctuations of <p according to the distribution 
((4]). The random stress tensor, which is directly related 
to the LB noise variables fj, is a Gaussian white noise 
source with correlations given by 



r/(r) ( 5 ai 5 f 3s + <W<Vy - -<5 Q( 3<5 7 <5 ) + C(r) S a p5js 



8(r - r')5(t - t') 



(74) 



Above expression is identical to the standard Landau- 
Lifshitz result for ideal fluids [9(|, except for the locally 
varying viscosities required to take into account possi- 
ble spatial inhomogeneities of the fluid. Note that any 
additional error terms in the Navier-Stokes eqs. (|69| orig- 
inating from the LB model have been neglected. These 
terms generally depend by a positive power on the den- 
sity gradient or flow velocity |79f and are expected to be 
negligible in the present case. 



B. Setup 

Simulations in the critical regime require fine-tuning of 
parameters as implied by the phase-diagram (Fig Q} as 
well as by LB-specific constraints. First of all, since the 
flow velocity must not exceed the lattice sound speed a s 
(where a 2 s = 1/3 for D2Q9 models), expression (pj?]) for 
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the fluctuation variance, k B T/po = 
k B T < <T 2 p ■ 



directly implies 
(75) 



Second, the density must remain strictly positive. Due 
to the Gaussian character of the density fluctuations this 
implies that the average density fluctuation should re- 
main much smaller than the mean density. Neglecting for 
simplicity spatial correlations, density fluctuations have 
a variance of (Ap 2 ) = p$k B T / c 2 ; hence, requiring that 
(Ap 2 ) 1 / 2 < po leads to 



k B T < c s p . 



(76) 



which is a more stringent constraint than eq. (fT5"j). since 
c 2 ~ x ^ a s f° r a non-ideal fluid in the critical regime. 
Thus the fluctuation temperature T must be chosen suf- 
ficiently small for the velocity fluctuations not to vio- 
late the approximate incomprcssibility of the LB method. 
Typically, values of k B T = 10" 7 1 .u. or less are sufficient. 

In the critical regime, it is particularly important to en- 
sure accurate equilibration of the fluid at all scales, since 
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FIG. 4: Equilibration of momentum in the critical region. 
j x denotes the i-component of the momentum, while and 
j± denote the projections of j longitudinal and transversal to 
the wavevector k. Simulation parameters (in l.u.): k, = 10~ 4 , 
r = 4.96 x 10~ 5 , u = 5.16 x 10~ 2 , k B T = 10~ 7 . All LB 
relaxation times are set to r = 0.8. 



here mode-coupling effects are dominant and thus errors 
induced at the smallest scales can propagate to larger 
ones and possibly infect the whole simulation. Since the 
noise covariance (|66p was derived in the limit of k — > 0, 
it is expected that using a sufficiently small value for 
k ensures equilibration to high accuracy. In fact, it is 
found that values of k < 10 -3 are already sufficient in 
the present case. This is demonstrated in Fig. 01 where 
- for a set of parameters in the critical region (see be- 
low) - the Fourier-transform of the spatial correlation 
function of the fluid momentum j = pu is compared to 
the theoretically expected result, (|j Q (k)| 2 ) = p k B T. To 
achieve a high statistical accuracy, correlations have been 
computed by averaging over 5000 snapshots over a total 
simulation time of 10 7 timesteps. As the figure shows, 
perfect equilibration at all scales with an error below 2% 
is obtained. 

The chosen value of k leads to restrictions on the pos- 
sible values of the free energy parameters r and u in the 
critical region, as these quantities enter the reduced f and 
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u [eq. (fT2j)]. which define the phase-diagram (Fig. Q} of 
the Ginzburg-Landau model. In order to observe Ising- 
type critical behavior, one would want to avoid too close 
proximity to the displacive limit and hence choose a large 
value for f [l|,[26||. However, as outlined in 9l|, for the de- 
terministic LB method there exists, as a consequence of 
the finite-difference approximations to the spatial deriva- 
tives, a lower limit for the interface width of around 1 l.u., 
below which LB simulations produce potentially wrong 
dynamics even for small density differences. Since the 
large-scale fluctuations at the critical point essentially 
consist of phase-separated regions ("critical droplets" ) 
that continuously break-apart and reconnect [92j, it is 
plausible that the restriction on the interface width de- 
rived for the non-fluctuating case continues to hold in 
some form also for fluctuating critical domains. Since the 
interface width of the critical domains is roughly given 
by the mean-field result 



y/-2n/r = 



(77) 



one obtains an upper bound on f of around 2. This ar- 
gument implies that one can not get arbitrarily close to 
the Ising limit without sacrifying correct dynamics. Sim- 
ulation results obtained in this work, however, indicate 
that this appears to be not too severe of a restriction 
for the application of the LB method to critical phenom- 
ena. Taking, for instance, ksT = 10~ 7 , w ~ 0(1) and 
K ~ O(10 -4 ), a simulation in the critical regime then re- 
quires that f c ~ 0(1), u c ~ O(l), implying r c ~ O(10~ 4 ) 
and u c ~ O(10 _1 ). To obtain a more precise location, 
one can gradually change one of the parameters r, u, 
or T and visually inspect the density field, compute the 
structure factor or apply a finite-size scaling analysis, as 
shown below. Note that, in principle, one can traverse 
the critical regime by either changing T, u or r, keeping in 
each case the other parameters fixed. To stay in line with 
the usual field-theoretical notion of the tcmpcraturc-likc 
variable, only r will be varied in the present work. 

In the critical region, relaxation processes become ex- 
tremely slow (critical slowing down), requiring, in par- 
ticular, at long wavelengths, a large simulation time in 
order to collect a sufficiently large number of statisti- 
cally independent samples (cf. [HI). At criticality, sound 
waves are overdamped (to be discussed in more detail in 
Ull) and decay with a rate of T(fc) = c 2 (k)/vi, where 
the generalized speed of sound is given by (see, e.g., 
[lU) c 2 s (k) = c 2 + pnk 2 and vi = (i] + 0/P i s the lon- 
gitudinal viscosity for a 2D fluid. Consequently, the 
largest relaxation time of the order-parameter can be 
estimated as t p ~ vijc\ (fc m i n ) ~ S 2 vi/(c 2 fi + 4ir 2 pK,), 
with S being the system size, fc m ; n = 2n/S the mini- 
mum wavenumber and c 2 denoting the value of c 2 for 
a correlation length of £ = 1. Here, the critical (mean- 
field) scaling of the thermodynamic speed sound sound, 
c 2 cx 1/x = c 2 £~ 2 oc 5~ 2 , has been used. Thus, the per- 
formance of a simulation can be optimized by choosing a 
small value of v\. (Note that this, however, also leads to 
an increase of the momentum relaxation time, which is 



cx S 2 /v.) For instance, S = 256, k = 1CT 3 and v x = 1CT 2 
gives a density relaxation time of t p ~ 10 6 l.u. Hence, 
a simulation must run roughly 10 s timesteps until accu- 
rate statistical information (errors less than one percent) 
for order-parameter related quantities is obtained when 
averaging over a few hundred realizations. This require- 
ment of simulation time might seem excessive, but one 
should keep in mind that, due to the global conservation 
of the order-parameter, large- wavelength fluctuations re- 
quire the rearrangement of mass over large distances. 



IV. RESULTS 

A. Correlation function 

In a simulation, the structure factor is computed from 
the order-parameter field <p on the discrete lattice by 



O(k) = — 

V ; N 



E^ r ) gik ' r 



(78) 



where the brackets indicate time-average over many sta- 
tistically independent samples and N is the total num- 
ber of lattice points. Due to periodic boundary condi- 
tions and the real-valuedness of </), it suffices to consider 
the structure factor in the first half of the first Brillouin 
zone, i.e. in the wavenumber range < k < n. Global 
mass conservation enforces 0(0) = (this point is ex- 
cluded from the plots). Fig. [S] shows the structure factor 
together with sample snapshots of the order-parameter 
field above and at the critical point. 

Above the critical point (Fig.[S£i), the correlation func- 
tion assumes a simple Ornstein-Zernike form [cq. (1251) ]. 



0(k) = 



1 + k 2 i 2 



(79) 



with £ being the correlation length and \ the compress- 
ibility. Note that in the above equation, k 2 re pre sents the 
Fourier-transformed discrete Laplacian (see [16[), which 
reveals itself in a deviation of the high-fc-part of C(k) 
from a simple k~ 2 power-law expected in the continuum 
case [l6j]. The discrete lattice effect becomes noticeable 
for wavenumbers k > 1. Note that around r = 0, self- 
energy corrections suppress the correlation length and 
compressibility below their mean-field values given by 
cqs. (J2SJ), (j2T)|) . Treating £ and x as n t parameters in 
cq. (|79|) allows one to extract the critical growth of the 
correlation length and compressibility upon approaching 
r c from above (see below). 

At the critical point, the correlation function is ex- 
pected to assume a power-law, C(k) ~ k~ 2+v [eq. ([T7|) |. 
for k > l/£. In the above equation, r) is the anomalous- 
dimension exponent, which takes a value of r\ = 1/4 for 
the 2D-Ising universality class. In Fig. [SJj, the structure 
factor close to criticality is shown for two different sys- 
tem sizes of 96 2 and 512 2 lattice sites and the same set 
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FIG. 5: Correlation function C(k) vs. wavenumber k obtained from a simulation (a) far above and (b) close to the critical point. 
The insets show the representative order-parameter field. In (a), the dashed curve represents a fit to an Ornstein-Zernike form 
[eq. (|25|) ]. while in (b), it represents the critical power law fc _2+ ' 7 . In (b), the structure factor for two different system sizes, 
.5=256 (•) and 96 (□). Note the finite-size effects at low k. Simulation parameters: (a) r = 10 - ,u = 4x 1(T 3 , k = 1CT 4 , 
k B T = 10" 7 ; (b) r = -10 -4 , « = 4,k = 2.0x 10~ 5 , k B T = 5.31 x 10" 10 . 



of simulation parameters. It is seen that while the criti- 
cal power-law behavior is well obtained in both cases in 
the intermediate wavenumber range, the structure fac- 
tor shows quite pronounced finite-size effects at small k. 
Specifically, for the larger system, the Ornstein-Zernike 
type 'shoulder' at low k indicates that the system is still 
slightly above its critical point, whereas the \ow-k excess 
seen for the smaller system apparently suggests that the 
system is already sub-critical. However, according to the 
results of a finite-size scaling analysis performed on the 
order-parameter and susceptibility (see below), not only 
the lager, but also the smaller system is still above its 
critical point. Although finite-size effects appear to be 
more pronounced in the present case as compared to, for 
instance, molecular dynamics simulations [12j . it is well 
known that different quantities (e.g., structure factor and 
order-parameter distribution) in general show a differ- 
ent finite-size scaling behavior and are governed by their 
own apparent critical points 0, HtJ ■ Here, for all stud- 
ied parameter combinations and system sizes it is found 
that the thermodynamic quantities are critical when the 
structure factor already displays slight effects of phase- 
separation (i.e., has an excess at low k). It is clear from 
Fig. [5}j that these discrepancies can in principle be re- 
duced by using larger systems. However, the convergence 
appears to be quite slow in the present case. Simulations 
performed deeper in the Ising regime furthermore sug- 
gest that this behavior is not directly associated with 
the proximity to the displacive limit. It might, however, 
be possible to improve the simulations by determining 
an optimized set of coupling constants in the free energy 
functional (so-called "perfect action^) for which scaling 
corrections are most suppressed [l], [9J, IH| . 

The fact that sufficiently far above the critical point 
the structure factor assumes an Ornstein-Zernike form 
allows one to extract the critical growth of the correla- 
tion length £ and compressibility \ by fitting expression 
([75)1 to the simulation data for C(k). As Fig. [5] shows, 



the correlation length and the compressibility approach 
the critical point by power-laws £ ~ 6~ u and x ~ ^ 1 
with exponents that asymptotically agree with 2D-Ising 
values v = 1 and 7 = 7/4. Further away from the criti- 
cal point (6 > 0.2), we observe cross-over to mean-ficld- 
like power-law behavior. It should be remarked that, 
especially in the case of the correlation length, the data 
admits in fact a certain range of values for the expo- 
nent v and critical temperature r c . For instance, in the 
present case it is found that the correlation length can 
be equally well described by an exponent of v rj 0.8 and 
a slightly different r c . Similar "effective exponents" have 
also been observed in previous Monte-Carlo simulations 
of the 4 -model [26| and reflect the fact that the width of 
the asymptotic region, where Ising-type behavior is ob- 
served, depends on the proximity to the displacive limit 
pj. Indeed, approaching a critical point that is located 
closer to the "Ising-limit" on the critical line is found to 
already restrict the possible fit values for v to a narrower 
margin around 1. Due to the finite size of simulation 
box, however, it is not possible to follow the correlation 
length up to arbitrarily small reduced temperatures 6. 

In the crossover regime from mean-field to critical be- 
havior, it is interesting to compare the simulation re- 
sults with the predictions of perturbation theory (section 
III A5P . Fig. [7J shows the correlation length as extracted 
from Ornstein-Zernike fits to the structure factor ver- 
sus the inverse of the dimcnsionless coupling constant 
A -1 = rn/uT, varying here only r. We see that, for 
A -1 ^> 1, non-linear effects are negligible and the cor- 
relation length closely follows the mean-field prediction 
£ = (w/r) 1 / 2 (dashed line). Once A becomes of the order 
of unity, fluctuation corrections to mean-field behavior 
grow, leading to a suppression of the correlation length 
from its mean-field value (which diverges at A -1 = 0). 
The solid curve in Fig. [7J represents the prediction for 
the renormalized correlation length £ = (ft//) 1 / 2 ob- 
tained from the numerical solution of the self-consistent 
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FIG. 6: Supercritical growth of (a) the correlation length £ and (b) the susceptibility in dependence of the reduced temperature 
6. £ and \ are extracted from fitting Ornstein-Zernike functions [eq. (|79|l ] to the correlation function C(k) [eq. J75J] obtained 
from simulations. Solid lines represent power-laws with 2D Ising exponents, while dashed lines are power-laws with mean-field 
exponents. 
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FIG. 7: Dependence of the correlation length on the dimen- 
sionless coupling constant A -1 = m/uksT in the crossover 
regime from mean-field (A 3> 1) to critical behavior (A ~ — 1). 
The dashed curve represents the mean-field correlation length, 
while the solid curve shows the prediction of perturbation 
theory, obtained from the numerical solution of eq. (|34[> . 
The symbols represent the correlation length extracted from 
Ornstein-Zernike fits to the structure factor obtained from 
simulations. 



equation (|34|) for the parameter r' . We see that the simu- 
lation results for £ agree well with the theoretical predic- 
tions until the onset of critical growth (at A -1 « —0.6). 
The eventual breakdown of the simple perturbation the- 
ory close to the critical point is not surprising since all 
wavenumber-dependent contributions - which are strong 
in 2D - to the renormalized parameters are neglected. 
To increase the accuracy of the predictions, more sophis- 
ticated renormalization group methods would have to be 
employed [95| . 



B. Finite-size scaling 

Before turning to the finite-size scaling analysis at 
the critical point, the general dependence of the basic 
thermodynamic quantities on temperature and system- 
size in the wider critical region shall be investigated. 
Fig. [5] shows the evolution of the block order-parameter 



Ml, susceptibility xl and specific heat cjj.l f° r different 
coarse-graining lengths L in dependence of the tempera- 
ture, which is here represented by the inverse dimension- 
less coupling A^ 1 = rn/uT (k, u and T are fixed). The 
curves are drawn as a guide to the eye [they represent 
the finite-size scaling functions defined by eq. (|50l) ]. 

Passing from the high- to the low-temperature phase 
(i.e. decreasing A -1 ), the order-parameter (Fig. [5^) dis- 
plays a sudden increase at around A -1 ~ —1.65, which 
can be identified with the critical point. Due to the def- 
inition of Ml as the average of the absolute value of the 
coarse-grained order-parameter in each subbox [eq. (|45[l ]. 
Ml is non-zero even in the disordered regime, but ap- 
proaches zero with increasing subbox dimension L. Al- 
ternatively to eq. (|45|) . the order-parameter can be de- 
fined by the position of the maximum of the underly- 
ing distribution Pl (inset to Fig. [5^,) , in which case the 
order-parameter is exactly zero in the disordered phase 
(except in the immediate neighborhood of the critical 
point, where, in 2D, the order-parameter distribution 
develops a bimodal structure, cf. sec. IIV C[) . In the low- 
temperature phase, fluctuations decrease the value of the 
order-parameter with increasing coarse-graining length. 
Due to inevitable interfacial contributions in the coexis- 
tence regime, this effect is more pronounced for Ml de- 
fined through eq. (|45p . For intermediate coarse-graining 
lengths L, the susceptibility - which for simplicity is com- 
puted in Fig. |5p via the same eq. (|T7|) both in the high- 
and low-temperature phase - shows a peak at the appar- 
ent critical point, consistent with the behavior of Ml- 
For very small L, the peak is "smeared out" to a shoul- 
der, while for the largest L (not shown in the plot), the 
low-temperature data is strongly affected by interfacial 
contributions. Note that the peak positions are practi- 
cally independent of L. The specific heat (Fig. [Sfc) shows 
a rapid increase around the critical point A ~ —1.65, but 
no peak, in contrast to Monte Carlo simulations [U [28| . 
The behavior of the specific heat seems to be similar 
to the order-parameter Ml and the absence of a peak 
might thus be related to the presence of interfacial con- 
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FIG. 8: Temperature dependence of (a) the coarse-grained order-parameter Ml [eq. (|45|) ]. (b) the susceptibility \l [eq. (|46fl ] 
and (c) the specific heat [eq. (|49[) ] in the critical region. The inset in (a) shows the position of the peak of the order-parameter 
distribution Pl (for visibility, only three curves are shown). The temperature is represented by the inverse dimensionless 
coupling A -1 oc r. The curves in each plot correspond to different subsystem dimensions L: • 1, ■ 2, ♦ 4, A 8, T 16, o 32. Data 
for L > 5/8 has been excluded. In all cases, u = 2.8 x 1(T 2 , k = 9.6 x 1CT 5 , k B T = 1(T 7 . 



tributions in the sub-critical regime. It should be finally 
remarked here that the order-parameter, susceptibility 
and specific heat obviously depend systematically on the 
coarse-graining length L. Hence, in order to obtain their 
true values, they have to be extrapolated to L — > oo. 
This does not affect the finite-size scaling behavior and 
is discussed further in sec. I IV CI 

In the immediate vicinity of the critical point, the the- 
oretical correlation length exceeds the system size and 
one has to perform a finite-size scaling analysis to ex- 
tract critical properties. First, the scaling with subbox- 
size L of the block order-parameter, susceptibility and 
specific heat at the critical point is investigated, keeping 
the total size S of the simulation box (here, S = 256) 
and all other simulation parameters fixed. In this case, 
the FSS ansatz (|5T))) predicts for the order parameter Ml 
and susceptibility xl the scaling forms 

M L ~ L-^»g M (L/S) , 
XL ~ L^g x (L/S) , (80) 
ch.l ~ log(L)g c (L/S) , 

where gj\/, g x and g c are scaling functions with lim- 
its go(L/S) — s- for L ps S due to the global order- 
parameter conservation and go(L/ S) — > const, for L suf- 
ficiently smaller than S. Note that the temperature de- 
pendence has dropped out of the above scaling forms, 
as the temperature is kept fixed at its presumed critical 
value. As Fig. |H] shows, for L < S/8 the simulation re- 
sults for Ml and Xl agree well with the FSS predictions 
of eq. (|80|) for the 2D-Ising case (solid lines in the plot). 
In case of the specific heat, the expected logarithmic scal- 
ing is obtained for L > 4 and extends to block sizes up 
toLsi 5/4. It is remarked that this close agreement can 
be only obtained in a rather narrow range around the 
critical point. Scaling plots like Fig. fallow in principle 
to estimate the true value of various intensive quantities 
by simple extrapolation of the straight line fits to the full 
system size L — S. 

Figure [TUJ shows the finite-size scaling behavior of 
the coarse-grained order-parameter and susceptibility for 



varying subbox sizes L and temperatures 8. Data for 
L = 1 as well as L > S/8 have been excluded from the 
plots, as they are expected to lie outside of the regime 
of validity of the FSS ansatz eq. ([50)1 . In the plots, 2D 
Ising exponents are used for the scaling transformation, 
together with a value for the critical temperature that 
is identical to the one obtained from the previous analy- 
sis. The supercritical (6 > 0) branch of the susceptibility 
(filled symbols in Fig. [TUb) shows good scaling behavior, 
while the scaling collapse is less satisfactory for the super- 
critical branch of the order- parameter (Fig.fTOb). For the 
latter quantity, this behavior persists for all tested values 
of the critical temperature r c . Note that the definition 
of the coarse-grained order-parameter [eq. f|45[) ] implies a 
non-zero value for Ml even in the disordered regime. Ml 
slowly approaches zero with larger temperature 9 or sub- 
system size L. For the sub-critical branches (open sym- 
bols), neither the data for the order-parameter nor the 
susceptibility collapse onto a master curve. Several rea- 
sons for the apparent scaling violations might be in order: 
First, the -model is equivalent to the Ising model only 
asymptotically close to the critical point and the width of 
the asymptotic region gets smaller with decreasing dis- 
tance to the displacive limit (r — > 0, u — >• 0) fil [26|. 
In fact, it is well-known that the FSS form ([50)) repre- 
sents only the leading order term of the full FSS expres- 
sion P, [HH, with the leading corrcction-to-scaling term 
being given by L~ u go{L l / v 6) (uj = 4/3 in 2D). Thus, 
corrections to scaling are necessarily always present in a 
simulation. If a higher level of accuracy is desired, one 
might perform simulations using an optimized set of cou- 
pling constants in the free energy functional, for which 
the leading-order scaling correction due to the dominant 
irrelevant operator is absent [H,|93|, 94]. However, the gain 
in using an improved set of parameters might be spoiled 
by the presence of the additional scaling variable L/S, 
which in turn requires a rather large size S of the total 
simulation box. Second, due to the global conservation of 
the order-parameter, domains of equal amounts of liquid 
and vapor are coexisting below the critical point, lead- 
ing to pronounced interfacial contributions to the order- 
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FIG. 9: Finite-size scaling of the block order parameter Ml (a), susceptibility %l (b) and specific heat ch,l (c) with the 
block size L (at fixed S) at the critical point. The solid lines represent the FSS predictions for the 2D-Ising case. For 
comparison, scaling laws according to mean-field case are shown by the dashed lines in (a) and (b). Simulation parameters are 
r = -4.8 x 1CT 5 , n = 2.8x 10~ 2 , k = 9.6 x 10~ 5 , k B T = 10" 7 . 
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FIG. 10: Finite-size scaling plots of the subbox order-parameter Ml (a) and susceptibility \l (b). The exponents /3, 7 and 
v are fixed to 2D-Ising values. The critical temperature r c defining is set to r c = —4.82 x 10" °. Legend: * 8 — 0.042, i 
6 = 0.0028, ♦ 6 = 0.013, ■ 9 = 0.006, • 6> = -0.001, « = -0.016, □ 9 = -0.044. Data for L = 1 and L > 5/8 have been 
excluded. 



parameter distribution, which deteriorate scaling in the 
ordered phase. Similar effects have been pointed out in 
the context of lattice gas simulations [33| 6 . 



Finally, the usefulness of the FSS ansatz written in the 
form (|5ip to locate the critical point is demonstrated. In 
Fig. 1 1 1 1 the appropriately rescaled order-parameter and 
susceptibility versus the non-linear coupling u is plot- 
ted, keeping all other simulation parameters fixed. By 
eq. (|5ip . the intersection point of all curves can be iden- 
tified with the critical point = 0, which, for the present 
choice of simulation parameters, occurs for a value of 
u ~ 6.8 — 7.0 x 10~ 3 . As expected, this value slightly 
depends on the quantity under consideration, but is oth- 
erwise consistent with the estimates of the critical point 
location from the FSS analysis of Fig. |9] 



In standard grand-canonical Monte Carlo simulations, interfacial 
effects are much reduced as one stays in a pure phase most of 
the time |53|| . 



C. Order-parameter distribution 

In the previous section, the FSS behavior of averaged 
thermodynamic quantities at the critical point was inves- 
tigated primarily. We shall now turn to a more detailed 
study of the behavior of the order-parameter and suscep- 
tibility in the off-critical regime as well as the properties 
of the underlying probability distribution. 

Far above the critical point (Fig. [L2k). the order- 
parameter distributions have a perfectly Gaussian shape 
centered around the mean order-parameter value 0. The 
variance decreases from the smallest block size (L = 1) 
towards the largest L = 5/2, which is understandable 
from the fact that coarse graining the system over a scale 
L integrates out fluctuations on smaller scales, which 
then do not contribute anymore to the variance. As dis- 
cussed in section III B[ the coarse-grained susceptibility 
XL, as determined by the width of Pl, in general dif- 
fers from the true susceptibility obtained in the thermo- 
dynamic limit due to the neglect of correlations at the 
boundary of the subsystem. In particular, in the off- 
critical regime (£ < L), \l is expected to differ from the 
true susceptibility by a correction factor ~ 1/L. This is 
demonstrated in Fig. 112b . where \l is computed from 
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FIG. 11: Determination of the critical coupling by finite-size scaling analysis according to eq. (JSTJ. (a) Plot of the subsystem 
order-parameter Ml [eq. Q45|)] scaled by L®' u vs. the coupling u. (b) Plot of the subsystem compressibility xl [eq. (|46p ] scaled 
by L - ' 1 '" vs. the coupling it. In all cases, r = —4.8 x 10 -5 and ksT = 10~ 7 and 2D-Ising values of the exponents /3, 7 and 
v are used. Different symbols correspond to subsystem sizes of L — 2 (o), 4, (a), 8 (□), 16 (o), 32 (•). Data for L = 1 and 
L > S/8 are excluded from the plot. The curves are drawn as a guide to the eye. 
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FIG. 12: Order-parameter distribution Pl (a) and corresponding reduced susceptibility Xl/x O 3 ) f° r different L far above the 
critical point (the correlation length is £ = 1). The individual curves in (a) correspond to different coarse-graining lengths L, 
where L — 1 for the outer curve and L — S/2 for the innermost curve. For better visibility, individual data points are not 
shown. The subbox susceptibility xl is extracted from Pl via Gaussian fits (•), from the variance of Pl (□), and from the 
central height -Pl(O) (a). The expected susceptibility x is computed from perturbation theory. The dashed line in (b) has a 
slope of -1. 



Pl by three different methods. As expected, the rela- 
tion xl ~ 1/i holds in the range ( < L « S, while for 
L ~ S, xl bends down towards zero due to the global 
order-parameter conservation. By extrapolating the lin- 
ear part in l/L towards L — > 00, the true susceptibility 
can be estimated. Good agreement between the extrapo- 
lated value and the susceptibility x obtained from mean- 
field theory is observed. 

In the critical regime (£ > L), the corrections due to 
missing boundary correlations will clearly not be given 
anymore by a simple surface-to-volume ratio as above. 
Instead, the susceptibility will gradually approach its 
FSS form xl ~ L 1 ^, as can be seen in Fig. Q21 where 
Xl obtained from the variance of Pl slightly above the 
critical point is plotted against L (cf. Fig. [3]). In the fig- 
ure, susceptibility data for L > S/8 is neglected in order 
to prevent a possible spurious influence due to the large 
value of L/S, for which the curves start to bend towards 
zero. Plotted in this way, one finds that extrapolating 
Xl to L — > S agrees well with the susceptibility obtained 
from Ornstcin-Zcrnikc fits to the structure factor of the 



entire system (Fig. [5]). 

Distinctly below the critical point, the order-parameter 
distribution is characterized by two displaced Gaussians 
centered around the spontaneous order-parameter value 
(Fig. ITla). Note that the probability distribu- 
tion covers more than three orders of magnitude between 
its center and its peak. The width of each Gaussian peak 
decreases with larger coarse-graining length L as more 
and more fluctuations are averaged out. The region be- 
tween the peaks arises from interfacial configurations and 
is significantly in excess of a pure Gaussian contribution. 
In agreement with the heuristic arguments outlined in 
section IIV CI the central region does systematically in- 
crease with larger subsystem size L and becomes approx- 
imately flat, as is expected by the presence of two-phase 
configurations with arbitrary proportions of liquid and 
vapor in each subbox. However, even for the largest sub- 
box sizes, the peaks are still far more dominant than the 
central region. This can be explained by two facts: first, 
interfacial free energies are still not negligible compared 
to bulk contributions, which, however, would be required 
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FIG. 13: System size dependence of the coarse-grained suscep- 
tibility xl obtained from the variance of the order-parameter 
distribution Pl slightly above the critical point (where L < 
£ < S). The dashed lines are extrapolations to the to- 
tal system size L — S and thus to the true susceptibility. 
The crosses represent the susceptibility obtained from the 
Ornstein-Zernike fits to the structure factor (cf. Fig. [(J). At 
the critical point, the susceptibility is expected to scale with 
the subbox size as L 7// " (dotted line) (cf. Fig. [SJJa). Different 
symbols correspond to reduced temperatures of 9 = 0.854 (•), 
0.270 (□), 0.124 (♦), 0.0219 (o). Data points for L > S/8 are 
excluded from the plot. The lines are drawn as a guide to the 
eye. 



in the thermodynamic limit. Second, and more impor- 
tant, for deep quenches, the liquid region (which in the 
present case is a single extended stripe) is not moving ap- 
preciably during the simulation time and thus each sub- 
box will be mostly covered by the same, virtually static, 
phase configuration. This is also indicated by the strong 
irregularities found in the central region of the distribu- 
tions for large L. To obtain the correct coarse-grained 
distribution, one would additionally have to perform an 
average over different simulation runs, which is, however, 
not attempted here. 

In principle, the thermodynamic order-parameter Ml 
can be defined either by the average over half of the distri- 
bution, Ml — (|tol|) [cq. (l4"5"1) ]. or by the position m max 
of the maximum of Pi(m). In the coexistence regime, 
one finds that the latter definition is in general in closer 
agreement with the theoretical prediction, eq. (|38p . for 
all values of L (Fig. ri4"b). In principle, a slight sys- 
tem size dependence is always expected as fluctuations 
in general tend to reduce the average order-parameter 
value, which is clearly seen closer to criticality (inset 
to Fig. [5Ji and Fig. \Wk). The order-parameter defined 
by eq. (|4"5j) strongly decreases with larger L due to con- 
tributions from phase-separated states to the average of 
\rtiL |. Thus, far above or below the critical point, defining 
the order-parameter Ml as the location of the maximum 
of Pl seems in general preferable over the definition of 
eq. (|45]) , since the former ensures that Ml is exactly zero 
in the high-temperature phase and has a negligible de- 
pendence on system size or interfacial contributions in 
the low-temperature phase. In contrast, the definition of 



cq. (|45[) behaves smoother in the critical region and is 
therefore better suited for finite-size scaling analyses. 

Analogously to the high-temperature case, one can 
obtain the coarse-grained susceptibility either from the 
peak height [\l ^ L d /87rTP|(m max ), see eq. {57)] or 
the peak variance. As all these methods implicitly as- 
sume the presence of two displaced Gaussians [cq. {S7J], 
the central region of the distribution should be excluded 
beforehand 7 . Due to the significant asymmetry of the 
wings and the pronounced interfacial contributions, the 
variance is thus most reliably obtained by fitting Gaus- 
sians to the peaks. It is seen from Fig. [Fib , that all 
the three different estimates of the susceptibility roughly 
agree, except for values of L close to the total system 
size. Extrapolating the linear part in 1/L of the coarse- 
grained susceptibilities \l to the limit L — > oo allows to 
obtain the true susceptibility. 

In Fig. [15] the order-parameter and susceptibility in 
dependence of the temperature (represented by the in- 
verse dimensionless coupling A -1 oc r) are compared to 
the predictions of perturbation theory (see sec. Ill A 5|) . 
The order-parameter data shown in Fig. [T5b is obtained 
from the location of the peak of the distribution, which 
is roughly independent of L and has negligible statistical 
scatter (see Fig. 114b) . One sees that the simulation re- 
sults for the order-parameter agree well with the predic- 
tion of eq. ([55|) for (<fi) including leading-order fluctuation 
corrections. As expected, deviations become noticeable 
closer to the critical point (here, A" 1 « —1.5), where 
a perturbative treatment is not applicable. One further 
notes that fluctuations generally lead to a decrease of the 
order-parameter value over its mean-field value, even far 
away from the critical point. The susceptibility data \ 
shown in Fig. [T5b is obtained by extrapolating the block 
susceptibility xl obtained from the peak height of the 
distribution to L — > oo (cf . Fig. (Hb) . In contrast to the 
order-parameter, the susceptibility data exhibits signifi- 
cantly stronger statistical scatter, in particular, closer to 
criticality. Nevertheless, for the temperature range in- 
vestigated, acceptable agreement between simulation re- 
sults and the prediction of perturbation theory for x i n 
the symmetry broken phase, eq. (|41[) . is found. 

In Fig. 1161 the coarse-grained order-parameter distri- 
bution at the critical point is shown. Interestingly, for 
sufficiently large coarse-graining lengths, the distribu- 
tion shows a pronounced double-peak structure, which 
is found to persist even slightly above the critical point. 
This is in agreement with Monte- Carlo results of the two- 
dimensional Ising model [52| and renormalization group 



calculations [71]. For small coarse-graining lengths, 
where the distribution essentially probes non-universal 
properties, Pl depends significantly on the location on 
the critical line: close to the displacive limit, P\ appears 



7 In particular, the normalization should be computed with the 
central region set to zero 
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(a) "I ( b ) LIS ( c ) 1/L 

FIG. 14: Order-parameter distribution Pl (a), order-parameter Ml (b) and susceptibility xl ( c ) for different coarse- graining 
lengths L far below the critical point. The individual curves in (a) correspond to different L, where L = 1 for the outer curve 
and L — S/2 for the innermost curve. In (b), the order-parameter as defined by eq. (|45[1 (•) and by the peak position of Pl (■) 
is shown. The dashed line represents the expected value of Ml from perturbation theory, eq. (|38[) . In (c), the susceptibility \l 
is extracted from Pl via Gaussian fits (•), from the variance of Pl (□), and from the height of the peaks PL(m max ) (a). The 
true susceptibility x is computed from perturbation theory, eq. (|41[) , The dashed line in (b) has a slope of -1. 




FIG. 15: Temperature dependence of the thermodynamic order parameter M (a) and susceptibility \ (b) in the phase- 
coexistence regime. The temperature is expressed here in terms of the inverse dimensionless coupling A -1 cx r. The symbols 
(•) represent simulation results, the solid curve represents the predictions of leading-order perturbation theory [eqs. (|38[) and 
(|41l) ] and the dashed curve the prediction of mean-field theory. The critical point is located here at A -1 ~ —1.5. 



concave, whereas towards to Ising-limit, it develops a 
double-peak structure. This is understandable since the 
gradient term in the free energy functional dominates 
over the bare Landau potential for high degrees of dis- 
placiveness. The scaling ansatz ([59| for the critical order- 
parameter distribution predicts that when expressing the 
data in terms of the scaled variables mL^ v and PlL - ' 3 /", 
all points should collapse on a single curve. In order to 
compare our results with the Ising model calculations of 
ref. [52j j . the rescaling procedure is implemented here by 
multiplying our data appropriately by the standard de- 
viations (to 2 ) 1 / 2 , which is expected to be equivalent con- 
cerning the overall scaling behavior since (to 2 ) ~ L -2 ' 9 /". 
It was checked that both ways to perform the rescaling 
lead to similar results. As Fig. H~6b shows, the predicted 
scaling of the distribution function expressed by eq. (|59p 
holds in a range 1 <C L <C S. This might seem surprising 
insofar, as the finite-size scaling of the low-order moments 
of P L (the coarse-grained order-parameter and suscepti- 
bility) works well already for the smallest box sizes (sec 



Fig. [5]). However, the full probability distribution obvi- 
ously contains more information than just its low-order 
moments, and thus a scaling of P L is only a sufficient, 
but not necessary condition for the scaling of Ml and 
Xl- br fact, in the present case, the fourth-order cumu- 
lant already fails to show the well-known scaling behav- 
ior observed in standard (grand-canonical) Monte-Carlo 
simulations of the Ising- model [52| . A similar behavior 
has also been observed in lattice gas simulations with 
a conserved order-parameter [33| . The scaling behav- 
ior of the distribution for smaller coarse-graining lengths 
is found to improve with increasing distance from the 
displacivc limit. A direct comparison of Pl in the scal- 
ing regime to the corresponding order-parameter distri- 
bution of the two-dimensional Ising model at criticality 
[53 | . represented by the solid points in Fig. 116b 8 , shows 
close agreement, except for a slight underestimation of 



8 To match the width of the distributions used here, the Ising 
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FIG. 16: Order-parameter distribution Pt(m) for different L at the critical point. In (a) the raw data is shown, while in (b) the 
distributions are rescaled by their standard deviation to compare with the predictions of eq. H59[) . For better visibility, curves 
in (b) that deviate from the scaling prediction are shown dashed. The filled circles in (b) represent the scaled distribution 
function from simulations of the 2D Ising model at the critical point by [52] . The curves in each figure correspond to different 
coarse-graining lengths L, where in (a), L = 1 for the outer curve and L = S/2 for the innermost curve; vice versa in (b). In 
all cases S = 256. Note also the linear axes scale, in contrast to Figs. 1121 and 1141 



the peak heights. 

V. SUMMARY 

In this work, static critical phenomena of a one- 
component fluid have been studied using a fluctuating 
non-ideal gas LB model. In this model, the fluctuating 
Navier-Stokes equations for the density and momentum 
of a compressible, isothermal fluid based on a Ginzburg- 
Landau </> 4 -free energy functional are solved. It is found 
that the model is able to capture the essential features 
of the static critical behavior associated with the 2D 
Ising universality class. A characteristic property of the 
present simulation method is the global conservation of 
the order-parameter, which demands a more careful in- 
terpretation of finite-size scaling results compared to, 
for instance, standard grand-canonical Monte-Carlo ap- 
proaches. The conserved nature of the order-parameter 
leads to the presence of coexisting two-phase states be- 
low the critical point and is expected to be the main 
source of scaling corrections in the present case. Despite 
these complications, the critical behavior of the struc- 
ture factor, order-parameter, susceptibility and specific 
heat is found to be overall well reproduced. However, it 
was noted that finite-size effects appear to have a quite 
strong effect on the structure factor, which assumes the 
expected critical scaling law at a slightly higher temper- 
ature than the other thermodynamic observables. For 
future work, it would be interesting to compare these re- 
sults to other LB models of non-ideal fluids. The order- 
parameter distribution function, which shows scaling be- 
havior at criticality and contains useful information on 



two-phase states below the critical point, compares well 
with theoretical predictions and Ising model calculations. 
Also, issues relevant to coarse-graining and generic fluc- 
tuation induced effects on observable quantities near and 
far from the critical point have been discussed. 

The present work has only dealt with static critical 
phenomena. While it is clear that Monte Carlo methods 
are usually better suited for this task, an assessment of 
the LB method in this regard is nevertheless important 
since the successful reproduction of equilibrium aspects 
is a necessary prerequisite for a faithful application of the 
method to, for instance, dynamical problems. Also, un- 
derstanding the equilibrium behavior of the model and its 
coarse-graining properties are important for many prac- 
tical problems employing an effective free energy descrip- 
tion, such as nuclcation and spinodal decomposition. The 
present work is thus hoped to provide a useful starting 
point for further applications of the LB method to prob- 
lems of current interest involving phase-transitions and 
critical phenomena in fluids. 
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